Spin- and charge-density waves in the Hartree-Fock ground state of the 

two-dimensional Hubbard model 
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The ground states of the two-dimensional repulsive Hubbard model are studied within the unre- 
stricted Hartree-Fock (UHF) theory. Magnetic and charge properties are determined by systematic, 
large-scale, exact numerical calculations, and quantified as a function of electron doping h. In the 
solution of the self-consistent UHF equations, multiple initial configurations and simulated anneal- 
ing are used to facilitate convergence to the global minimum. New approaches are employed to 
minimize finite-size effects in order to reach the thermodynamic limit. At low to moderate inter- 
acting strengths and low doping, the UHF ground state is a linear spin-density wave (1-SDW), with 
antiferromagnetic order and a modulating wave. The wavelength of the modulating wave is 2/h. 
Corresponding charge order exists but is substantially weaker than the spin order, hence holes are 
mobile. As the interaction is increased, the 1-SDW states evolves into several different phases, with 
the holes eventually becoming localized. A simple pairing model is presented with analytic calcu- 
lations for low interaction strength and small doping, to help understand the numerical results and 
provide a physical picture for the properties of the SDW ground state. By comparison with recent 
many-body calculations, it is shown that, for intermediate interactions, the UHF solution provides 
a good description of the magnetic correlations in the true ground state of the Hubbard model. 

PACS numbers: 75.30.Fv, 71.15.Ap, 71.45.Lr, 71.10.Fd, 75.10. Lp, 75.50.Ee 



I. INTRODUCTION 

The Hubbard model is one of the most fundamental 
models in quantum physics. Despite numerous analytic 
and numerical investigationsir— , key questions still re- 
main about the properties of this model^i"— . Surpris- 
ingly, even at the mean-field level, its phase diagram has 
not yet been fully determined, and the ground state mag- 
netic properties arc not completely known. 

The Hubbard model was originally proposed to de- 
scribe correlations between d-electrons in transition 
metalsi^. At half- filling (one electron per lattice site), it 
gives a simple description of the so-called Mott insulator, 
with antiferromagnetic (AFM) order. Soon after the dis- 
covery of high- Tc cuprate superconductors, it was pointed 
out that the two-dimensional (2D) Hubbard model might 
be an appropriate minimal model for high-Tc cupratesi^, 
because of the copper-oxygen plane geometry and the 
proximity of the superconducting transition to the AFM 
phase of undoped mother compounds. The 2D Hubbard 
model has since become a focal point of research in con- 
densed matter and quantum many-body physics. 

Recently, rapid experimental progress in optical lat- 
tice emulatorsi^ has promised a new way of to approach 
Hubbard-like models. Using ultra-cold fermionic atoms 
trapped in periodic laser-field potentials, these highly 
controllable experiments are capable of potentially 'sim- 
ulating' the Hubbard model directly. Thus the properties 
of the Hubbard model are not only of importance theo- 
retically but can also be of direct experimental relevance. 
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The Hartree-Fock (HF) method is the simplest 
paradigm to describe a quantum many-fermion system. 
The method finds the single Slater determinant wave- 
function which minimizes the variational energy. As is 
well known, the mean-field approximation involved can 
turn out to be very severe. Nevertheless the HF method 
has often provided the foundation for our qualitative un- 
derstanding of many systems in condensed matter and 
quantum chemistry. For example, HF correctly predicts 
an AFM order in the ground state of the Hubbard model 
at half-filling, even though the strength of the AFM order 
is overestimated and translational symmetry is (neces- 
sarily) broken. In quantum chemistry, HF is the starting 
point for most calculations and serves as the basis for 
understanding the electronic structure of many systems. 

Because correlation effects (e.g., the correlation energy, 
which is a fundamental concept in the framework of den- 
sity functional theory )i^ are often defined using the HF 
solution as a reference, qualitative and quantitative un- 
derstanding of the HF state is of key importance. This 
has not always been easy to achieve. For example, the 
nature of the unrestricted HF (UHF) state in the electron 
gas at high and intermediate densities was only recently 
detcrminedi^. 

In the Hubbard model, HF calculations are in princi- 
ple straightforward. The 2D Hubbard model has been 
studied within the HF approximation in some of the 
pioneering works on high-Tc superconductors. Inhomo- 
gencous states have been found at small dopings, such 
as spin polaronsii, domain walls^^^— , and spin den- 
sity waves^i"— (SDW), and phase diagrams have been 
proposedii^i^i. Due to computing power limitations, 
however, these studies have either done exact numer- 
ical calculations at only a few doping and interaction 
parametersiSi^SiS^, or have scanned parameters with re- 



strictcd forms of the solutioni^i^^i^. Furthermore, finite- 
size effects were difficult to remove, as we discuss below, 
which can mask the true solution in the thermodynamic 
limit. A systematic and quantitative understanding of 
the magnetic properties of the UHF ground state has 
not been achieved. 

In this work, we perform extensive numerical calcu- 
lations to determine the exact UHF ground state of 
the Hubbard model in the low to intermediate inter- 
acting strength regime. The exact UHF ground state 
we achieved is a full numerical solution of HF Hubbard 
Hamiltonian (see Sec. [H] Eqs. ^ and (|3])), as opposed 
to constrained searches or non-self-consistent solutions. 
We study the spin and charge properties as a function of 
interacting strength and doping concentration. Full nu- 
merical solutions of the UHF equations are computed us- 
ing twist-averaged boundary conditions for system sizes 
well beyond those previously studied. We also present a 
simple pairing model, with analytic calculations at low 
doping and small interacting strengths, to complement 
the numerical results and provide a qualitative physical 
picture of the magnetic properties of the model. 

Our combined numerical and analytical calculations 
show that, at a finite doping h, the UHF ground state 
at low and intermediate strengths U/t is a static lin- 
ear SDW (1-SDW) state. As the interaction strength 
is raised beyond a critical value, 1-SDW order develops 
along the [10]-direction, accompanied by a weaker linear 
charge density wave (1-CDW). The characteristic wave- 
length of the 1-SDW is found to be 2/h and the wave- 
length of the corresponding 1-CDW is 1/ft,. As the in- 
teraction strength is increased, stripe or domain walls 
states develop along the diagonal [ll]-direction, in which 
the holes are localized. The diagonal stripe (d-stripes) 
state and the 1-SDW state are separated by either a lin- 
ear stripe state (1-stripes, along [10]-direction) or a di- 
agonal SDW (d-SDW) state, depending on the doping. 
These are summarized with a UHF phase diagram for 
interaction up to U/t ^ 10 and doping up to h ^ 35% 

The remainder of the paper is organized as follows. In 
Sec. [ni the self-consistent scheme used for solving the 
mean-field Hubbard model is summarized. The numeri- 
cal results are presented in Sec. IIIIl and analytic calcula- 
tions are described in Sec. IIVI In Sec. |V] the results are 
discussed and summarized in a phase diagram, and we 
conclude the paper in Sec. IVII 



II. METHOD 

The Hamiltonian of the single-band repulsive Hubbard 
model reads 

"H = -t ^ [cl^Cr'a + cl^Cra) + C/ ^ '^rt'^r;, (1) 
{rr'},(T r 

where U > is the interacting strength and t is the hop- 
ping amplitude between nearest neighbor sites (denoted 



by {rr'} in the summation). Throughout this work, en- 
ergy is quoted in units of t and we set i = 1. The opera- 
tor cl^ (cra) creates (annihilates) an electron with spin a 
{a =t, i) at site index r, which runs through the lattice of 
size N = Lx X Ly. The total number of spin-cr electrons 
is denoted by N„, and we assume that the system has no 
spin polarization, i.e. N^ = N^. Under this assumption, 
the model has only two parameters, namely, the onsite 
repulsion U and the doping 



l-(Ar^ + iVj/iV = iVhoic/iV, 



(2) 



where we have used A^hoio to denote the number of holes 
in the system. Due to particle-hole symmetry, we confine 
ourselves in the region where N^ + N^ < N. Therefore 
the total density is given by (n) = 1 — h. 

Standard linearization of Eq. ([l} leads to the mean- 
field HF Hamiltonian: 



nuF = n 
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with 



^&F = -* E (4.Cr'. + cl^C, 



{rr'} 



+ UY,nra{nra)-^UY,{nrt){nrl), (4) 



where a is the conjugate of a and (rira) is an average 
density. The mean-field decoupling employed in Eq. ([3]) 
assumes the z-axis as the quantization direction, thus 
breaking the spin rotational symmetry of the Hubbard 
Hamiltonian ([T]). After fixing the quantization orienta- 
tion and requiring no spin polarization, the solution of 
the HF Hamiltonian is restricted to the S^ = sector, i.e. 
spin textures in the x-y plan, for instance spiral SDWs, 
are excluded. (At low U, the solutions turn out to be 1- 
SDWs. Then a single spiral cannot be the ground state, 
since a left-handed spiral can always be combined with a 
right-handed one, or vice versa, to make an 1-SDW which 
has lower energyJ^) 

For a given set of parameters {U,N,N^,Ni), the 
HF Hamiltonian (|3]) is numerically solved using a self- 
consistent scheme. We begin the procedure by selecting 
a trial solution in the form of a single Slater determinant 
for each spin component: 
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where each column is normalized to 1. In the restricted 
HF (RHF) method, spin-f and spin- J, parts of the to- 



tal wavefunction are the same: $. 
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The RHF 



method always gives the non-interacting solution in the 
systems studied in this paper. In the UHF method, which 
is adopted here, $i and $, ' are allowed to differ and 
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they converge via the coupled Eqs. (|4]). The trial densi- 
ties at site r can be expressed as 
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where 'i?' indicates conjugate transpose of the matrix, 
and we have assumed that the orbitals in $Jj are or- 
thonormal. An A^ x iV matrix M^ (M^) for nl^p CH^f) 
is then constructed from the densities. By exactly diag- 
onalizing Af^, we obtain the energy 



Ei' 



Vi 









4=1 



^(1) 



(7) 



< A^^ are the lowest N^ 
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eigenvalues of M'^. The wavefunction ^a is obtained 
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by filling up N„ corresponding orbitals of A^ 
density (nU 

which is used to update M^ (-M^)- We iterate this pro- 
cess until the total energy E^^^ = E. 
density {rirJ 



W 



e[''' and the 

is converged. 

Care must be taken when updating the density during 
the iteration. As is typical in self-consistent algorithms, 

convergence to a fixed point is not guaranteed if (nj-o- ) 
is taken directly as an input for the £-th step. To improve 
convergence, we adopt a mixing scheme: The i-th input 
density is constructed as a linear combination of previous 
input and output densities as: 
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where 'in' indicates the input density to construct M'^, 
and 'out' denotes the output density calculated by diago- 
nalizing Af^ . The mixing parameter a is typically chosen 
to be between ^ 0.5 and 0.75. 

Due to non-linearity of the coupled Eqs. (|3]), we imple- 
ment two additional procedures to help the system reach 
the global minimum. Firstly, different initial wavefunc- 
tions arc used and the consistency between the results is 
checked. Secondly, we perform multiple annealing cycles: 
in each cycle a random perturbation (whose strength can 
be controlled) is applied to the converged solution and 
the self-consistent process is repeated. 

To reduce shell and one-body finite-size effects, we use 
twist-averaged boundary conditions (TABC)^"— , under 
which the wavefunction \I'(ri, r2, . . .) gains a phase when 
electrons hop around lattice boundaries 
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where L is the unit vector along L, and the twist angle 
= i^xj^y) is an input parameter which is randomly 
chosen in this work. For a given 0, the TABC is the 
same as a random shift of the momentum space grid. 
This reduces the discretization error in the integration. 
In the HF solution, the TABC is applied to each orbital. 
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FIG. 1: (Color online) Contour plots of CD (left) and SD 
(right) at half-filling and U = 4.0 on a 16 x 16 lattice. CD 
is uniformly distributed at a density of 1. SD is AFM with 
uniform amplitude. 



i.e., each column in Eq. ^. With a generic 0, there 
will be no degeneracy in the one-electron energy levels. 
We often average the results over many random twist 
angles^! in each system to improve convergence to the 
thermodynamic limit. As can be seen from the energy re- 
sults in Sec. IIIIl this procedure produces a smooth curve 
vs. doping, where the one-body finite-size effect is inin- 
imized. Additional finite-size errors, which result from 
the interaction and the formation of long wavelength col- 
lective modes2^, are not removed from this approach. We 
use rectangular lattices in our simulations to help detect 
the 1-SDW states with long modulating wavelengths, as 
discussed in Sec. IIIIl 



III. NUMERICAL RESULTS 

Various observables are computed with the converged 
UHF wavefunction. Two quantities examined through- 
out the paper are the charge-density (CD) p{r) and the 
spin-density (SD) s(r) defined as 



p(r) = (71rt) + (nr^), 
S(r) = (jlrf) - {Url). 



(10) 
(11) 



We will also study the converged UHF eigenvalues Ao-p 
and momentum distribution ripo- ~ {c^ Cp^r), where Cpo- 
is, as usual, defined by the Fourier transform of Cra- Fig- 
ure [T] illustrates the behavior of the reference system at 
h = 0, which is an AFM state with constant p{r) = 1. 

In Fig. [21 CD in the 16 x 16 reference system is plotted 
as holes are doped into the lattice. As doping h is varied, 
holes tend to cluster and form different patterns. These 
patterns have a strong ft,- dependence, which is a result 
of strong finite-size effects. Here the system is at an 
intermediate interaction strength of f7 = 4.0. As the 
interaction becomes weaker, we find that the variations 
in the patterns become larger and depend sensitively on 
(not shown). This is similar to what is seen in the 
UHF solution of an electron gasi^. 
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FIG. 2: (Color online) Contour plots of CD for a supercell of 
16 X 16 at [/ = 4.0 as Nhoie is increased. The finite-size effect 
is strong until the 1-CDW wavelength is decreased sufficiently 
to fit into the simulation cell, as shown in the right bottom 
plot. 



FIG. 4: (Color online) Contour plots of CD for systems of 
fixed Af = La; X Ly = 32 X 32 = 16 X 64 = 8 X 128 = 1024 (from 
top to bottom) and fixed doping h = 3/32 {Nhoie = 96) at 
U = 3.0. A stable 1-SDW solution emerges when the supercell 
is commensurate. Note that only the accompanying CDW is 
shown here. 
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FIG. 3: (Color online) Contour plots of SD for a supercell 
of 32 X 32 at [/ = 3.0 as NhoU is increased. An 1-SDW exists 
when the density is such that the supercell size is sufficient 
to accommodate the 1-SDW, as shown in the left top or the 
right bottom plot. 



A. Linear spin-density wave (1-SDW) state 

We first focus on low to moderate interacting strengths 
[U < half of the bandwidth) and small doping {h < 0.1), 
and examine the properties of the UHF solution as a 
function of doping h, i.e., as the system moves away from 
half- filling (ft, = 0). We will show that the UHF ground 



state at low and moderate C/ is a linear spin-density wave 
(1-SDW) along the [01] direction. Figure [3] shows the re- 
sults from a 32 X 32 supercell. An 1-SDW is seen whenever 
the density is such that an 1-SDW can be accommodated 
in the supercell. (The choice between x- and y-dircctions 
in the broken-symmetry UHF state is of course random. 
To help visualization in the figures, we have selected the 
same direction, either by an initial bias or by rotating the 
final result.) At incommensurate densities, strong finite- 
size effects are present, where the pattern of the cluster is 
not scalable to the thermodynamic limit. An example is 
seen by comparing iVhoio = 16 in Fig. [5] (not long enough 
for one period of SDW) and A^hoie = 64 in Fig. |31 in both 
cases h = 1/16. iVhoie = 24 in Fig. [5] vs. iVhoic ~ 96 in 
Fig. |3]is another (both have h = 3/32). The finite-size 
effects will be further discussed below. 

Although significantly larger lattice sizes are reached, 
the pattern variation clearly indicates that care must be 
taken in a numerical calculation, and additional ingre- 
dients arc needed, in order to better approach the ther- 
modynamic limit. Wc use two additional ingredients in 
our numerical simulations: TABC and rectangular super- 
cells. To reduce the one-body finite-size effects, most of 
our results are averaged over ~ 20 random values. In 
plots showing 0-averaged results, the statistical uncer- 
tainties from the twist angles are indicated by the error 
bars. The residual (two-body) finite-size effects are re- 
duced by the use of rectangular supcrcells. This allows us 
to study longer wavelength modes without increasing the 
computational cost (compared to a square lattice of the 
same number of lattice sites, N.) Obviously rectangular 
lattices break the symmetry between x- and y-directions, 
and can introduce an additional bias. To minimize the 
effect, we carry out calculations with different supercells 
with varying aspect ratios to check consistency in the 
results. 

An illustrative set of results is shown in Fig. SI We 
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FIG. 5: (Color online) Contour plots of CD (left) and SD 
(right) vs. doping. The system is an 8 x 64 supercell at 
U = 2.0 with doping of /i = 1/32, 2/32 and 3/32 (from left to 
right). The wavelength of the 1-CDW is Ai_cdw = l/h and 
that of the 1-SDW is Ai_sdw = 2/h. 



X Ly — 



adjust Lx and Ly while keeping the size N = L^ 
32 X 32 fixed. An 1-SDW solution is seen in a rectangular 
supercell whenever Ly is sufficiently large to accommo- 
date a wave. Note that the rectangular supercell does 
not bias the SDW in the y-direction (when Ly > L^,). 
An 1-SDW is observed along the i-direction if Lx is com- 
mensurate with the SDW wavelength. (An example of 
this is in Fig. [6] below, where the solution in the 20 x 36 
lattice is two waves propagating along x-direction.) 

From the results in Figs. [2] and |3l it is clear that the 
wavelengths of the 1-SDW and 1-CDW vary with doping 
h. The results of an 8 x 64 lattice with various values of 
h are shown in Fig. [5] As can be seen, the wavelength 
of the 1-CDW/SDW decreases with h. Unlike in Fig. H 
the lattice size in this case has been chosen so that Ly 
is commensurate with the wavelength in each figure. For 
example, there are exactly two CD waves at h = 1/32, 
giving a wavelength of Ly/2 = 32(= 1/h). The wave- 
lengths of SDW (right panel) are twice those of CDW. 
When the doping is doubled or tripled, the number of 
waves being accommodated changes accordingly, i.e. the 
wavelength shortens by 1/2 or 1/3, respectively. The 
modulating wavelengths of the 1-CDW and 1-SDW are 
thus given by 



Al-CDw(/i) 
Al-SDw(/l) 



1 

h' 
2 

h' 



(12) 
(13) 



The wavelength relations are verified with many different 
choices of the aspect ratio. 

The variational energy of the UHF ground state is ex- 
amined in Fig. [B] A series of supercells are studied with 
a fixed N — LxX Ly — 720 and h = 0.1, while varying L^ 



FIG. 6: (Color online) Ground state energy per particle as 
a function of the aspect ratio for a series of supercells with 
a fixed N = L^ x Ly = 720. Doping is at /i = 0.1 and the 
interaction strength is [/ = 2.5. Results are averaged over 22 
random values; statistical error bars are shown, although 
some are too small to be seen. For the supercells which can 
accommodate full 1-SDW/CDW, whose wavelengths are de- 
termined by Eqs. (|12p . H13p . the variational energy is consis- 
tent and lower. 



and Ly. It is seen that, for all supercell choices commen- 
surate with the predicted wavelength, the energies are 
consistent and are lower. In systems which are incom- 
mensurate and cannot accommodate the 1-SDW/CDW, 
the resulting ground state energies from the UHF solu- 
tion are higher, indicating the frustration effect in the 
variational solution because of the finite size of the su- 
percells. In Sec. IIVI we will present an analysis showing 
why in general the 1-SDW is favored at low U . 

More lattice sizes at various dopings and interacting 
strengths are studied. The amplitudes of the 1-SDW in 
the obtained solution are summarized in Fig. [T] It can 
be seen that at each fixed density, the 1-SDW amplitude 
decreases as U is decreased and eventually vanishes, in- 
dicating the disappearance of the broken-symmetry UHF 
solution at a critical interaction strength Uc- Below Uc 
only a RHF solution exists. The critical value Uc appears 
to decrease with h and approaches at zero doping. This 
is consistent with the situation at half- filling {h = 0), 
where the Fermi surface (FS) is an open shell and a UHF 
state can be formed by 'pairing— across it with no cost 
to the kinetic energy. For a fixed U , the amplitude of the 
1-SDW decreases with doping (as does the wavelength). 

The amplitude fluctuation is the strongest near Uc, in- 
dicated by large statistical errors, and decreases as U is 
increased. This can be understood from the mechanism 
for the 1-SDW states in the UHF solution. The 1-SDW 
state is formed by 'pairing' or nesting of electrons near 
the FS^ (see also Sec. |1V|. At low [/, the UHF so- 
lution only contains a small number of cxcitation o^^i^^ 
to plane-wave states immediately beyond the FS. In a 
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FIG. 7: (Color online) 1-SDW amplitude as a function of U 
at various doping for several supercells. At each data point, 
the result is averaged over 22 random values and the error 
bar is the statistical error. From left to right, the doping is 
increased. At a fixed doping, different supercells give consis- 
tent results. The amplitude increases with U beyond Uc and 
converges to a stripe or domain walls state (see Sec.lV)). 



FIG. 8: (Color online) UHF eigenvalues A^p vs. momentum 
p (top) and corresponding momentum distribution (bottom) . 
Both quantities are plotted along symmetry lines in momen- 
tum space, as depicted in the inset. The system is a 16 x 48 
supercell with doping oi h — 1/24 for a series of U. In the 
top, the RHF (non-interacting) band-structure is also shown 
for comparison. 



finite-sized system, how well the desired pairing can be 
achieved depends sensitively on the particular topology 
of momentum space grid, and the results therefore show 
more fluctuation with respect to iV or 0. Thus the 1- 
SDW amplitude is small around Uc, and sensitive to the 
boundary conditions, giving relatively large statistical er- 
ror bars. At larger U, there are more excitations above 
the FS, and the plane- wave states necessary for pairing 
become available independent of 0, so less fluctuation is 
seen. 

The picture we described above is supported by the 
UHF band structure and momentum distribution shown 
in Fig. [5] In the figure, we plot the UHF eigenvalues A-|-p 
(shifted by the mean-field background U{nr^) {n^^) /2) for 
a series of U . Each A-t-p is identified with a wavevector 
p by the maximum plane-wave component in the corre- 
sponding wavefunction, i.e., according to the magnitude 
of |(p|(/)^p)|. The corresponding momentum distribution 
is also shown. Results are the same for a =t and \, and 
are only shown for spin-f electrons. We will omit the 
CT-index below unless it is necessary. At small U values 
[U < 1/4 of the bandwidth), the deviation of rip^ from 
the non-interacting (or RHF) result is not drastic. We 
see that, as U exceeds Uc, a gap opens up in the band 
structure. Only a small number of states, |p), near the 
FS participate in the formation of the broken-symmetry 
state. As U is increased, there are more excitations and 
more states becoming involved. In Sec. llVCi we discuss 
the mechanism in further detail, and show how it is de- 
scribed by a simple pairing model at low U . 

As seen from Fig. [3 once the interaction strength is 
above the immediate vicinity of Uc, the finite-size effect 
becomes minimal in the system sizes we have studied. 



The wavelength and amplitude of the 1-SDW (CDW) do 
not change with the supercell size. Larger supercells give 
essentially identical results with the SDW replicated to 
fill the (conuirensuratc) supercell. 



B. Diagonal spin-density wave (d-SDW), linear 
and diagonal stripe (1/d-stripes) states 

As the interaction strength U is further increased, the 
UHF ground state changes character. Figure [5] shows the 
CD and SD along the y-direction, along which the linear 
wave propagates. Above Uc, the amplitude of the 1-SDW 
(and CDW) grows with U . As U is further increased, 
the CD reaches 1 and starts saturating, creating deeper 
density valleys at the nodes of the 1-SDW. The maxi- 
mum and minimum of CD and the 1-SDW amplitude as 
a function of U are plotted in Fig.[TU]to further illustrate 
this. As discussed in Sec. IIII Al CD/SD orders are de- 
veloped beyond U ^ 1.5. The 1-SDW amplitude is much 
greater than that of the 1-CDW. The CD maximum sat- 
urates at 1 above U ~ 2.5, indicating the formation of a 
linear stripe (1-stripes) state. The stripe or domain wall 
states differ from the SDW state because of CD satura- 
tion, forming hole-free domains that separate regions in 
which the holes are localized. The SDW state, in con- 
trast, is a wave state in which the CD spatially oscillates 
but does not reach 1, and the holes are dclocalized. 

Thus at low dopings (high densities, h < 0.1), the 1- 
SDW state turns into an 1-stripes state as U is increased, 
with the 1-stripes along the same direction {x- or y-) and 
having the same characteristic wavelength. When U is 
further increased the solution changes orientation, turn- 
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FIG. 9: (Color online) CD (top) and SD (bottom) along y- 
direction vs. U. The system being studied is an 8x64 supercell 
with doping of 1/32 at [/ = 1.0, 1.3, 1.5,2.0,4.0. Each curve 
is a ID cut in which the linear wave propagates. Beyond 
Uc, the 1-CDW and 1-SDW amplitudes increase with U and 
the ground state ends up in an 1-stripes state. The CDW 
amplitude is much weaker than that of the SDW. 
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FIG. 10: (Color online) Maximum and minimum of the CD 
and SDW amplitude for 8 x 64 supercell with doping of 1/16. 



ing into a stripes state with modulation along the [11]- 
direction, a diagonal stripes (d-stripes) state. 

At somewhat larger doping (0.1 < h < 0.3), the 
evolution of the 1-SDW state with U is different. The 
SDW state changes its modulation direction from the 
[10]-direction to diagonal. (d-SDW has been discussed 
in Ref. [22, for example.) Figure [11] shows an example 
for doping of /i = 1/6. We see that the modulating wave 
changes direction from [10] at f7 = 4.0 to [11] at [/ = 5.0, 
and the d-SDW saturates to become a d-stripes state at 
U = 9.0. 

We have scanned different parameter combinations to 
map out the sequence of the evolution of the UHF ground 
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FIG. 11: (Color online) Contour plots of CD (top) and SD 
(bottom) vs. interacting strengths. The system being studied 
is a 36 X 36 supercell with doping of fe = 1/6 at {/ = 4.0, 5.0 
and 9.0 (from left to right), representing 1-SDW, d-SDW and 
d-stripes state respectively. 



state. In Sec. [Vl a phase diagram is sketched to summa- 
rize the properties of the UHF ground state in the part 
of the phase space on which we have focused. The differ- 
ence in the pairing mechanism of the d-SDW state from 
that of the 1-SDW state is briefly discussed in Sec. IIVDI 



IV. ANALYTIC CALCULATIONS 

In this section we present a phenomenological model 
of the 1-SDW state at low U and small h. The model 
will help explain the numerical findings and provide a 
simple physical picture that captures the basic features 
of the exact UHF solutions in this parameter regime. 
The numerical studies are independent of the analysis 
here, but together they will give a more complete de- 
scription of the UHF states. Below we first discuss the 
basic pairing niodeU^— , then carry out calculations in 
detail in the limit of small U and h for the 1-SDW state, 
which is the focus of the present work. Some quantita- 
tive comparisons and validations of the pairing analysis 
are then presented, using the numerical data from calcu- 
lations presented in Sec. IIII Al We then briefly discuss 
the mechanism for d-SDW and d-stripes orders at higher 
U. 



A. Pairing model 

At low U , the region of interest in momentum space 
is the immediate vicinity of the FS, where pairing effects 
of electrons determine the nature of the UHF solution. 
(Often the effect has been discussed in the context of 
nesting. We refer to the mechanism as pairing since, al- 
though nesting greatly facilitates pairing in the Hubbard 
model, it is not required for the pairing mechanism to be 
realized, as seen in the electron gasi^.) In the fully filled 



71 







-71 



A, 






/^ 


"A 






Aq 




> ./ ^/ y^^ 1 


, -• 




X •■ y/^ / yr-- \l 


1 N 




/ ,#: ;io, -^v 


Ww 


\ 




\ 

•••■ / 


\ \ ^^■■. 


/ . 


/ 




V 


/ 


\ \ / X 




\ ■■•••...-•••■/ 






\ / 






\ / 













-71 







71 



FIG. 12: (Color online) Illustration of the pairing model at 
small I] and h. The half-filling FS is the large diamond (red 
dashed). The non-interacting FS at low doping remains ap- 
proximately the shape of a diamond (blue solid). AFM order 
arises from qo, the pairing vector across the half- filling FS. 
The pairing vector is q across the doped FS. The difference 
between qo and q, Aq, determines the characteristic modu- 
lating wavelength of the 1-SDW. 



region inside the FS, the electron density is uniform. 



"-pcr(p) = "TLcr 



N 
4^' 



(14) 



We first specify the pairing mechanisnJ^ more explic- 
itly. Recall that the non-interacting energy for the state 

|P) is 



Ep = — 2(c0SPa; -t- COSPj^). 



(15) 



The plane-wave state is |p) = ^Le'P'', with r = (x,y), 
where x and y arc integer coordinates denoting lattice 
sites. Consider a pair of spin-'f and spin-^ electrons in a 
p state, where ep < ep, with ep the Fermi energy. In the 
pairing modeli^, this pair is made to partially occupy a 
p' state: 



t) 



"pip; 
wp|p) 



^'plp;, 

^'pIp'), 



(16) 



where ImdP 



1. This pair gives the following 



contribution to the local density in real space 
1 



n^(r,p) 
%(r,p) 



N 
1 

N 



(l + 2|?/pWp|cos[(p'-p).r]), (17) 
(l-2|?/p«p|cos[(p'-p).r]). (18) 



And an SDW state will result from the state in Eq. (|16p , 
with local spin 

s(r,p) = n^(r,p)-n^(r,p) 

4 

== ]^l"p^'pl ■cos[(p'-p)-r]. (19) 

The SDW state lowers the interaction energy contribu- 
tion of the pair compared to the non-interacting solution 
(i.e., the solution when Vp = 0) by the amount: 



A£v(p) = UJ2 

r 

u 



n^(r,p)n^(r,p)-— — 



E^'(^'P) 



(20) 



If we have multiple pairs each formed as in Eq. (|16p , the 
change in interaction energy follows the same relation: 



U 



^^v^-^E 



U 



E 



E^^'^'P) 



(21) 



where the sum over p is over all pairing plane waves (one 
of the four sides is illustrated by striped areas in Fig. lT^ . 
At half-filling, the shell at the Fermi level, i.e. on the 
border of the diamond, is open, with the number of de- 
generate p states equal to twice the number of spin-f (or 
spin-|) electrons that need to be accommodated. Pair- 
ing can be achieved by choosing p' — p = qo = (tTjTt), 
i.e. having electrons occupy two states in the open shell 
across the FS. This is perfect nesting and the SDW 
formed has perfect AFM order. Because pairing occurs 
in the open shell at the FS, the reduction in interaction 
energy from the SDW has no penalty, i.e. no increase in 
the kinetic energy. 



B. The linear spin-density wave state 

We next consider the case of low U, slightly doped 
(U <^ bandwidth), to help understand the mechanism of 
the 1-SDW state. As the FS shrinks with small doping, 
we assume that it remains approximately the shape of 
a diamond. The distance between the FS at half-filling 
and the doped FS is determined by 



d — hn. 

4 



(22) 



As the interaction is turned on, it can become advan- 
tageous for some of the electrons near the FS to be par- 
tially excited. Partially occupied states around the FS 
can then allow pairing across the FS, which causes a cor- 
relation between electrons of opposite spins to generate 
an SDW. The presence of the SDW will lower the inter- 
action energy. However, in this more general case there 
will also be an increase in the kinetic energy. When the 
lowering of the interaction energy surpasses the increase 



in kinetic energy, an overall lower energy state is found 
compared to the free-electron (or RHF) solution. 

We first determine the kinetic energy change. At low 
U, pairing occurs near the FS. Electrons from a small 
region immediately inside the FS are excited. As a crude 
modeU^i^, we assume that a fraction / of the electrons 
within a distance S of the FS are excited, as illustrated 
by the horizontally striped region in Fig. 1121 The excited 
electrons occupy the region (vertically striped) immedi- 
ately above the FS, also of thickness ~ S. We take Up 
and Vp in the pairing state in Eq. (|16p to be independent 
of p: Up = u and Vp ~ v. Thus the vertically striped area 
has uniform density, and / = |wp . An upper bound to 
the kinetic energy increase due to this process is easily 
estimated. It is, for each excited electron, given by: 



A£k(p) 



VpEp • Ap, 



(23) 



where Ap ~ 6 (1, l)\/2/2. The total kinetic energy in- 
crease is then 



A^K = 8/n 
N5^ 



8/ 



AsKip) dS 
1 + cos 



(24) 



where S in the integral is over the horizontally striped 
area inside the FS, and the factor of 8 accounts for the 4 
sides and 2 spin species. 

We now determine the interaction energy change, and 
show that the optimal SDW is along the y- or x-direction. 
From Eq. (j2ip we see that the maximum reduction is 
achieved by maximizing the quantity 



^v-E 



E ^'^^^^ ■ 



(25) 



where the sum over q is over all pairing states, with 
q = p' — p. This is realized if all the electron pairs 
line up their pairing vectors. There are two groups of 
pairing states, corresponding to the two diagonal direc- 
tions. Within each group, the optimal choice is for all 
pairs to have one common pairing wavevector q. Let us 
denote the pairing wavevectors along [11] and [—11] by 
q and q', respectively, and write: q = (tt, tt) — Aq and 
q' = (— 7r,7r) — Aq'. We then obtain: 

/y oc y^[cos(q ■ r) + cos(q' ■ r)]^ 

r 

= 7V + ^(cos[(Aq + Aq')-r] 

r 

-f-cos[(Aq- Aq')T]). (26) 

The maximum is achieved in Eq. ([^5)) when Aq = ±Aq'. 
This occurs when q and q' are such that the SDW modu- 
lation from the two groups of pairing states are the same, 
leading to a positive 'interference' between them. The di- 
rection of the modulating wavevector must be along [01] 
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FIG. 13: (Color online) Absolute values of kinetic energy 
gained and interaction energy lost in the pairing model. On 
the left the energies are plotted as a function of h for several 
values of U . On the right, Uc is plotted vs. h. 



(or [10]). The magnitude is given by 
|Aq| = 2V2d = hn, 



(27) 



as illustrated in Fig. [121 This leads to the following total 
reduction in interaction energy: 



ASv 



.^.HV.(.-i 



(28) 



Thus the lowest energy state is an 1-SDW with broken x- 
y symmetry, with the modulation along either the x- or 
the 2/-direction. The modulating wavelength is Ai_sdw = 
2/h, consistent with our numerical result. 

To reach an SDW state of lower energy than the non- 
interacting solution, the condition 

\AEy\ > \AEk\ (29) 

must be satisfied. From Eqs. p4)) and |28)) . we obtain 



\u\'U 1- 



> 2 



1 + cos ■ 



hn 



(30) 



Taking \u\ ~ 1 on the left-hand side, we obtain a rough 
estimate to the critical value which U must exceed: 



C/n = 



hn 



h 

1 

2 



(31) 



The absolute value of the kinetic and interaction en- 
ergy changes in Eqs. ([24]) and (f28| are plotted vs. h in 
the left panel in Fig.[T31 AE'k is independent of U, while 
AEy is proportional to U, for which several curves are 
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FIG. 15: (Color online) Density plots of An(p), the momen- 
tum distribution difference from RHF solution (left) and its 
correlation An(p)An(p') (right). The system is a 16 x 24 
suppercell with doping of 1/12 at [/ = 3.0. Negative peaks at 
(±7r, ±(7r — 7r/12)) in the correlation result from the pairing. 



FIG. 14: (Color online) Energy plot of the modification to the 
RHF band structure in the UHF solution, for a sequence of U 
values in the 1-SDW regime. Shown are the values A^p — fp 
vs. p, where ep is given in Eq. (|15|) . The system is a 16 x 48 
supercell with doping of /i = 1/24. 



plotted for various values of U . It is seen that a critical 
value of U exists for doped system (ft, 7^ 0). Above C/c, the 
two curves cross at a critical /ic, below which the broken- 
symmetry 1-SDW state exists. As U increases, the point 
of crossing, he, moves to the right. Equivalently, the crit- 
ical C/c decreases as doping is reduced. In the right panel 
the curve of C/c vs. h is plotted to illustrate this. 



C. Comparison with numerical results 

The simple model and analysis above capture most of 
the properties of the exact UHF ground state at low U 
and small h. It gives the correct 1-SDW modulating wave- 
length, and explains the existence of C/c and how it varies 
with doping. Because of the crudeness of the model, 
the values of C/c and other quantitative features are not 
very accurate compared to the exact numerical results in 
Sec. IIII Al Larger discrepancies can be expected further 
away from its domain of validity, namely small doping 
and modest interaction (although it also incorrectly pre- 
dicts C/c = 4 as /i -> 0) . 

The model considers only pairing of two electrons, so 
CDW is excluded. This is consistent with the numerical 
result that at low C/, CDW is much weaker than SDW 
order. The exact UHF solution will necessarily involve 
more electrons in the pairingi^, which will lead to a larger 
energy lowering | AE'vl (and thus lower C/c) and will result 
in CDW, as observed in the numerical results. 

Figure [T3] shows the modification to the RHF band 
structure in the UHF solution as a function of interaction 
strength. The difference between the UHF eigenvalue 
Acrp and the RHF spectrum Cp in Eq. (fT5|) is plotted for 
all momentum values p. As discussed in Fig. [51 the eigen- 



value Ag-p is identified with the momentum p with which 
the corresponding cigcnstatc has the maximum overlap. 
We see that, just above C/c, a small fraction of the states 
on the FS are involved in pairing, which creates a small 
energy lowering that leads to the UHF solution. The plot 
is for a single twist angle. In a finite system, the shift in 
momentum space from the twist creates a small asymme- 
try between each pair of surfaces diagonally across. At 
small C/ > C/c, this is reficctcd in the solution as an asym- 
metry in the gaps on the two surfaces. As U increases, 
excitation spans a wider region at the FS, and the gap 
structure from pairing becomes more pronounced. 

The momentum distribution from a numerical UHF 
solution at h = 1/12 is shown in Fig. [151 The left 
panel plots ?T.(p) minus the non-interacting value 7io(p): 
An(p) = ^(p) — '^o(p)- Electrons arc excited from 
the darker area to the lighter. The right panel shows 
the two-point correlation function from the left panel: 
An(p)An(p') vs. (p — p'). Negative peaks are seen at 
(±7r, ±(7r — 7r/12)) on the right, which result from the 
pairing between the negative just inside the FS (where 
electrons are excited from) and the positive immediately 
above the FS (where electrons are excited to) in the left 
panel. The position of the negative peaks indicates a 
pairing vector of Aq = (0, /itt), consistent with the pair- 
ing vector in the analytical model. 



D. Diagonal spin-density Avave states 

As mentioned before, diagonal modulations lose the in- 
terference between [11] and [—11], so a diagonal (or any 
orientation other than [10] and [01]) SDW is not the so- 
lution at small h and moderate U . This does not exclude 
it as a solution as we move away from this parameter 
regime, when the distortion to the FS becomes more se- 
vere. 

This situation happens when the doped FS is deformed 
sufficiently away from the half-filling shape of a diamond 
and the area of excitation becomes sufficiently large to 
reach the half-filling FS. The number of pairs that could 
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FIG. 16: (Color online) Illustration of the pairing scheme 
for d-SDW order. The left panel shows n{p) and the right 
panel An(p), the difference from the non-interacting solution. 
The momentum distribution is actual numerical data from a 
system of 36 x 36 with doping of 1/6 at C/ = 5.0. Electrons 
are excited from the FS across the {—it, 7r)-direction to the 
FS across the other direction (7r,7r), such that the FS along 
the latter reaches the half- filling FS. This enables two groups 
of pairings to maintain interference, with Aq = Aq', to lower 
the energy. 



participate in the 'interference' of the 1-SDW is decreased, 
because the FS no longer has the shape of a diamond. 
Eventually it becomes energetically more favorable to 
have the FS be longer in one diagonal direction than 
the other, i.e., to break the four- fold rotational symme- 
try. As illustrated in Fig. [HI it is then possible to create 
two different types of pairing states along the two diago- 
nal directions, such that they share a common modulat- 
ing wavevector along one diagonal direction: Aq = Aq'. 
The two groups of pairs will achieve interference, similar 
to the case of 1-SDW. As in Sec. IIVI the pairing vector is 
determined by h, giving Aq = {hn, hir). which gives rise 
to an SDW with modulating wave along [ll]-dircction, 
and of wavelength Ad-sDW = \/2/h. The corresponding 
wavelength for d-CDW is l/y/2h. This is consistent with 
the numerical results in Sec. IIIIBI 



V. DISCUSSION 

We can now place our 2D results in the context of an 
HF phase diagram for the Hubbard model. Our numer- 
ical calculations have focused on small and intermediate 
dopings (h from to ~ 0.3), and small to moderate in- 
teractions {U from to ~ 10), because of possible con- 
nections with the many-body ground state at moderate 
interacting strengths. The analytic calculations are for 
small h and low U, where our pairing model captures the 
physics in the HF framework. Our numerical results arc 
sufficiently detailed such that we could determine some 
phase boundaries as shown in Fig. [T71 We fitted the 
numerical locations for the phase transition or crossover 
using power functions, except for the AFM to FM tran- 
sition which was fitted by an exponential. Because of 
the limited number of data points and the finite reso- 
lution with which the transition was scanned, there are 
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FIG. 17: (Color online) Phase diagram of the ground state of 
the 2D Hubbard model from UHF. The phase boundaries are 
determined by fitting our numerical results, and are meant 
only as rough guidelines. Solid lines separate the antiferro- 
magnetic (AFM) phase from the paramagnetic (PM) phase 
and the ferromagnetic (FM) phase. Within the AFM phase, 
the different regions include: 1-SDW (SDW state with a lin- 
ear modulation along [10]-direction); 1-stripes (density satu- 
ration to 1, with linear modulation along the [10]-direction); 
d-SDW (SDW state with a modulating along the diagonal 
[ll]-direction); d-stripes (density saturation to 1, with diago- 
nal modulation). The black dotted line gives the theoretical 
estimate (Stoner criterion) for the transition from the RHF 
solution (PM) to FM. 



significant uncertainties in the fits, of several line widths 
or larger. The phase boundaries are thus only meant as 
rough guides. 

At half-filling, the UHF solution is an AFM state. 
Upon doping, there is a phase boundary Uc{h), shown 
as the blue line in Fig. [T71 below which is the PM phase. 
Above Uc{h) is an AFM region where a rich set of sub- 
regions exhibit different characters, including the 1-SDW 
states we have focused on in this work; we describe this 
region in further detail below. Above the AFM phase 
is an FM phase. Our numerical UHF calculations show 
that the FM state has lower energy above the green solid 
line. The RHF approach, naturally, predicts an earlier 
transition to FM. This is the theoretical phase boundary 
from Stoner criterion, and is shown as the black dotted 
line. Recall that we have excluded spiral SDWs. As we 
discussed, this is not the ground state at low U (sec also 
Refs. [l9ll23r ). However, at large U, spiral orders can be- 
come more favorable deep in the d-stripes region. 

Between the PM and FM phases is the AFM phase. In 
this region, at low and intermediate U, we see an 1-SDW 
state with a long wavelength modulation along the [10]- 
direction; a weaker CDW accompanies the SDW. Near 
half-filling, as U is increased the 1-SDW state evolves 
into a 1-stripes state which shares the same characteristic 
wavevector as the 1-SDW, but whose CD saturates to 1 
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in regions separated by 'stripes' anchored by the nodal 
positions defined by the SDW. The holes are localized 
in these stripes. This is consistent with the observation 
in Ref. [22 of SDW deforming into domain walls with in- 
creasing U . The transition from delocalized holes (such 
as the 1-SDW state) to localized holes is denoted by the 
red dashed line in Fig. 1171 As we move further away from 
half-filling, the 1-SDW at lower interaction changes its di- 
rection of modulation as U is increased. This forms a d- 
SDW state. The transition from a state with modulation 
along the [10]-direction to one with diagonal modulation 
is denoted by the cyan dot-dashed line. We see that the 
two dash hues cross each other. At low doping (/i < 0.1), 
the system reaches an 1-stripes state first before chang- 
ing the direction of modulation to a d-stripes state. At 
higher doping, the order is reversed. The 1-SDW first 
changes into a d-SDW state. As U is further increased, 
density saturation appears, and holes become localized 
in a d-stripes state. 

It is important to keep in mind that the results we have 
discussed and the phase diagram above are for HF theory. 
For strong interactions in particular, the HF results are 
expected to be severely biased and correlation effects can 
fundamentally change the nature of the many-body state. 
For example, the FM phase was shown not to exist at low 
density {h > 0.5) in the 3D Hubbard model^S. 

The present work was in part motivated by a recent 
quantum Monte Carlo (QMC) calculation's, which indi- 
cated that the ground state of the 2D Hubbard model has 
a long wavelength SDW collective mode. Upon doping, 
the AFM order at half-filling was found to evolve into an 
SDW state with a long wavelength modulation which has 
essentially a constant charge-charge correlation at low to 
intermediate interacting strengths. Given that the UHF 
solution is qualitatively correct at half-filling, it was nat- 
ural to ask to what extent the UHF solution contains any 
of these features upon doping. 

We see from the numerical results in this work that the 
UHF solution appears to qualitatively capture the basic 
features of the magnetic correlations in the ground state 
upon doping, as it does at half-filling. Of course the UHF 
solution gives a static modulated SDW, while the many- 
body ground state in the QMC preserves translational 
invariance and the SDW correlation is only seen in the 
correlation functions'^. This is similar to the situation 
at half-filling. 

In the UHF solution, the tendency for the holes to lo- 
calize is much overestimated. This was part of the reason 
to focus on low U in the present study. A CDW corre- 
lation almost always accompanies the SDW in the UHF 
solution, and holes appear to localize (leading to domain 
walls or stripes) at C/ '--^ 4. In contrast, holes remain de- 
localized (wave-like) in the many-body solution's, with 
essentially constant charge-charge correlation, until the 
strong interaction regime {U > 10). It is an interesting 
question whether diagonal order, which is present in the 
HF solution at larger C/, is present in the true many-body 
ground state. 



The UHF solution thus provides a useful starting point 
for understanding the magnetic and charge correlations 
in the ground state of the Hubbard model at intermedi- 
ate interactions. In addition, the ability to reliably deter- 
mine the true UHF ground state numerically could prove 
valuable in QMC calculations, which often require a trial 
wavefunction and where the qualititativc correctness of 
the trial wavefunction can make a significant difference. 
Although the physics in the UHF solution is sensitive 
to the particular many-body Hamiltonian, the basic ap- 
proach we have used and the basic ideas of the analytic 
calculations are general (see also Ref. [Ig) and can be 
expected to find applications in other many-fermion sys- 
tems. 



VI. CONCLUSION 

In summary, we have performed exact numerical calcu- 
lations for the UHF ground state of the Hubbard model 
systematically for a wide range of lattice sizes, initial con- 
ditions, doping and interaction strengths. Special care 
has been taken to reduce finite-size effects in order to 
obtain the solution at the thermodynamic limit. These 
results allow us to map out the magnetic phase diagram 
for regimes most relevant in modeling condensed matter 
systems. 

A broken-symmetry UHF solution exists above a criti- 
cal Uc; whose value increases with doping. Above Uc{h), 
the ground state is a static 1-SDW/CDW, with a modula- 
tion whose wavelength is inversely proportional to doping 
at small h. The amplitude of the SDW/CDW decreases 
with h and increases with U. At low U, the SDW am- 
plitude is much stronger than that of the accompany- 
ing CDW; and the holes are essentially delocalized. For 
larger U, the SDW and CDW amplitudes become more 
comparable. At small doping, the solution turns into the 
1-stripes state with the same characteristic modulating 
wavevector and holes localized at the nodal positions, 
before eventually entering the d-stripes state. At larger 
doping, the 1-SDW state first turns into the d-SDW state 
before eventually entering the d-stripes state at larger 
interactions. 

We have also presented an analytic theory to explain 
the mechanism for the formation of the SDW state. The 
model provides a conceptual understanding of the physics 
of SDW which can be applied in systems beyond the 2D 
Hubbard model. Comparison with recent QMC results 
shows that the UHF solution captures the magnetic cor- 
relations in the true many-body ground state at interme- 
diate interactions. 
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